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1 Introduction 

Let f(x) be a real function that is rational except for the square root of a cubic or 
quartic polynomial with at least one pair of conjugate complex zeros. Then / f(x)dx 
can be expressed in terms of standard elliptic integrals with complex variables, which 
are subsequently changed to real variables by using quadratic transformations fToH [|TT|1 . 
Since the transformations complicate the formulas, it is desirable to have algorithms 
for numerical computation of standard elliptic integrals with complex variables. Such 
integrals are met in other contexts also, for example in conformal mapping, and the 
complex variables might not occur in conjugate pairs. 

This paper contains algorithms for numerical computation of complete and incomplete 
elliptic integrals of all three kinds when the variables are complex (with some restrictions 
for integrals of the third kind). They are similar to algorithms published earlier || for 
real variables, but several improvements that apply to the real case have been made 
in the course of extending them to complex variables. The integrals computed are the 
symmetric integrals of the first and third kinds and two degenerate cases, of which one 
is an elementary function and the other is an elliptic integral of the second kind. Other 
integrals can be obtained from these by using the formulas and references in Section 4. 
The method of computation is to iterate the duplication theorem and then sum a power 
series up to terms of degree five; the error ultimately decreases by a factor of 4 6 = 4096 
with each duplication. Since this method is slowest when the integrals are complete, we 
add a faster algorithm for computing complete integrals of the first and second kinds, 
including Legendre's K(k) and E(k) for complex k, by the method of arithmetic and 
geometric means. 

The symmetric integral of the first kind is 



where the square root is taken real and positive if x, y, z are positive and varies con- 
tinuously when x, y, z become complex. The integral is well defined if x, y, z lie in the 
complex plane cut along the nonpositive real axis (henceforth called the "cut plane"), 
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(1) 



with the exception that at most one of x, y, z may be 0. The same requirements apply 
to the symmetric integral of the third kind, 

3 r°° 

Rj{x,y,z,p) = - [(t + x)(t + y)(t + z)}~ l l 2 (t + P y l dt, (2) 
2 Jo 

where p ^ and the Cauchy principal value is to be taken if p is real and negative. A 
degenerate case of Rf that embraces the inverse circular and inverse hyperbolic functions 
(see Section 4) is 

Rc(x, y) = R F (x, y,y) = - (t + x)- l ' 2 {t + y^dt . (3) 

2 Jo 

It is well defined if x lies in the cut plane or is and if y 0; the Cauchy principal value 
is to be taken if y is real and negative. Professor Luigi Gatteschi pointed out to me that 
Fubini, while still a student at Pisa, proposed the use of a function equivalent to 1/Rc 
in his first published paper ||12|| . A degenerate case of Rj that is an elliptic integral of 
the second kind is 

3 r°° 

R D (x, y, z) = Rj{x, y,z,z) = - [{t + x)(t + y)}- 1/2 {t + z)~ 3 / 2 dt , (4) 



2 Jo 

which is well defined under the same conditions as Rf except that z must not be 0. 

Because Rd is symmetric only in x and y, it is sometimes convenient to use a com- 
pletely symmetric integral of the second kind, 

I roo /IX V Z \ 

Rc(x,y,z) = - 4 l [( t + *)( t + ,)( t + i )]-V^_ + ^ + _j trft , (5 ) 

where any or all of x, y, z may be and those that are nonzero lie in the cut plane. If the 
closed convex hull of {x, y, z} lies in the union of and the cut plane, Rq is represented 
by a double integral that accounts for its usefulness in problems connected with ellipsoids, 

I p2tt pit 

Rc(x,y,z) = — / / (x sin 2 9 cos 2 (j) + y sin 2 9 sin 2 (ft + z cos 2 Q) 1 ' 2 sin 9 dO d(f) . (6) 
47r Jo Jo 

Except that at most one of x, y, z may be 0, Rp has the same representation with exponent 

— 1/2 instead of 1/2. With the same exception and with x,y,z permuted so that z ^ 0, 

Rg can be obtained from Rf and Rd by using the relation 

1 I xy 

2R G (x,y,z) = zR F (x,y,z) - -{x - z){y - z)R D (x,y,z) + J — . (7) 

o V Z 



Algorithms for direct computation of Rc and Rf with real variables by successive Landen 
transformations are given in Jl . In Section 2 of the present paper, one of those algorithms 
is extended to complex variables, but only in the complete case. For each of the functions 
defined above, the complete case is the case in which one of x, y, z is 0. 

The functions Rp and Rc are homogeneous of degree —1/2 in their variables, Rj 
and Ro of degree —3/2, and Rc of degree +1/2. Their homogeneity and symmetry 
replace a set of linear transformations of Legendre's integrals, simplify their quadratic 
transformations and other properties, and make it possible (see f7fl~ IPII ) to unify many 
of the formulas in customary integral tables. 

The earlier versions of the algorithms in Section 2 (except the last one) were published 
in P|, modified in |6], (4), (5)] to avoid underflows, and implemented by Fortran codes in 
several major software libraries, in the Supplements to [|7| and ||, and in [ 13| and [14 



§6.11]. Codes in C can be found in [O, §6.11]. The algorithms in the present paper are 
preferable to those in || in the following four respects (in addition to incorporating the 
modification in |J): 

(1) All complex values of the variables are admissible for which Rf,Re>, and Rc 
were defined above, with the exception that computing the Cauchy principal value of 
Rc{x,y) when y < requires the preliminary transformation ([H]). The variables of 
Rj are restricted by conditions that are sufficient but not necessary to keep the fourth 
variable from being transformed to by the duplication theorem. 

(2) The algorithms have been rearranged to reduce substantially the number of arith- 
metic operations. 

(3) A bound on the fractional error of the result can be specified directly. 

(4) Because Rc is used repeatedly in the algorithm for Rj, its computation has been 
speeded up a little by including terms up to degree seven (instead of five) in the truncated 
power series. 

Section 3 contains information for checking codes based on the algorithms of Section 
2, and Section 4 relates other integrals to the integrals used here. 
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2 Algorithms 

We shall summarize briefly the method of computation; details are given in || for the case 
of real variables, and a full discussion of the complex case with proof of error bounds will 
appear eventually in a book now in preparation. Validity of algorithms can be tested in 
the complex domain by the consistency checks in Section 3 and sometimes by comparison 
of the complete case with the last algorithm in the present Section. 

If we obtain x m+ \ , y m +i , Zm+i from x m , y m , z m as prescribed by fllT]) below, the 
duplication theorem implies 



Although x m , y m , z m converge to a common limit in the cut plane as m — > oo, differ- 
ences like x m — y m decrease by a factor of only 4 when m increases by 1. When the 
fractional differences between x m , y m , z m and their arithmetic average become smaller 
than an amount determined by the desired accuracy of the final result, Rp is expanded 
in a multiple series of powers of these fractional differences, denoted by X, Y, Z. Because 
of the symmetry of Rp the series can be rewritten in terms of the elementary symmetric 
functions E\ , E 2 , E 3 of X, Y, Z; and because E\ = X + Y + Z = the terms up to degree 
5 are very simple. The truncation error, being of degree 6, ultimately decreases by a fac- 
tor of 4 6 with each duplication. An estimate of truncation error is provided by ||, (A. 10)]. 

ALGORITHM FOR Rp . We suppose that at most one ofx,y,z is and those that 
are nonzero have phase less in magnitude than it. The function Rp(x, y, z) defined by ([][) 
is to be computed with relative error less in magnitude than r. (We assume r < 3 x lO -4 .^ 
Let x = x, y = y, z = z, and 



R F {x , 



m = 1, 2, 3, . . . . 



(8) 



-4, 



x + y + z 



Q = (Sr)- 1 / 6 



max{|A - x\, \A Q - y\, \A - z\} . 



(9) 



— 
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For m = 0, 1, 2, ... , define 




A 



-m+l 



A m + A 
4~ 



■in 



(10) 
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^m+l — 7 j Urn+1 ~ ^ > z m+l ~ ^ > (,-L-U 

where each square root has nonnegative real part. 

Compute A m for m = 0, 1, . . . , n , where 4~ n Q < \A n \ . Define 

X = ^, Y = ^y, Z = -X-Y, (12) 
4M n 4M n ' K J 

E 2 = XY-Z 2 , E 3 =XYZ. (13) 

Then 

R F (x, y, z) « A7 n 1/2 (l - + ^ 3 + \-E\ - ^£ 2 £ 3 ) (14) 
with relative error less in magnitude than r. 



The statement that the relative error is less in magnitude than r means that the 
true value of the function lies inside a circle in the complex plane with center at the 
computed value and radius equal to r times the distance of the center from the origin. 
We assume that r is large compared to the machine precision, so that roundoff error 
is negligible compared to the error produced by approximations made in the algorithm. 
These remarks apply to all the algorithms in this paper. 

We note the relations 

^-m — g j ^im — X m — — , A — 1 — , l^lOJ 

and similar relations obtained by permuting (x, X), (y, Y), (z, Z) . Although Aq can be 
in the complex case, A n ^ because Q > 0. (It can be shown that A m ^ if m > 1.) If 
x m i Dm i z m are real and nonnegative, the inequality of arithmetic and geometric means 
implies \ m < 3A m , whence A m+ i < A m . 

If y = z then Rp reduces to Rc- Because E 2 and E 3 are no longer independent, the 
series in (|l4f) becomes a series in one variable. By including terms up to degree seven 
(instead of five), we can usually save one duplication, which seems worthwhile because 
computation of Rj requires one computation of Rc in each cycle of iteration. 
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ALGORITHM FOR Rc ■ Let x and y be nonzero and have phase less in magni- 
tude than 7T, with the exception that x may be 0. The function Rc{x,y) defined by (|3|) is 
to be computed with relative error less in magnitude than r. (We assume r < 2 x 10~ A .) 
Let Xq = x , yo = y , and 

A = ^-, Q = (Sr-r^Vo - x\ . (16) 

For m = 0, 1, 2, . . . , define 

= 2\/ X m y/y m + Urn i = ^ > (17) 

Xrn+1 = ^ , Vm+1 = T , (18) 

where each square root has nonnegative real part. 

Compute A m for m — 0, 1, . . . , n , where 4r n Q < \A n \ , and define 

y-A 



(19) 



n 



Then 



with relative error less in magnitude than r. 



If the second variable of Rc is real and negative, the Cauchy principal value is 

1/2 

, x + y 

by [|1(| (4.8)] and (|73|) . This vanishes if x — 0. 



Rc(x,-y) = ] R c {x + y,y), y>0, (21) 



In the notation of (|2T)| ) to (|2lf), the duplication theorem for Rj is 

1 6 

Rj{x m ,y m ,z m ,p m ) = -i2j(x m+ i,2/ m+ i,z m+1 ,p m+ i) + — l + e m ). (22) 

The duplication theorem for Rc has been applied to [|, (5.1)] to allow a wider range of 
complex phase for the variables. Iteration of (|22"D yields 

n—l ^—m 

Rj(x ,y ,z , p ) = 4~" Rj(x n , y n , z n , p n ) + 6 ^ — — 1 + e m ) . (23) 

m=0 dm 
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The first term on the right side can be treated as a symmetric function of the five variables 
x n , y n , z n , p n , p n by writing (t + p)~ l in (|2|) as (t + p) _1 / 2 (t + p)^ 1 ^ 2 . When the frac- 
tional differences X, Y, Z, P, P between these five variables and their arithmetic average 
become small enough, Rj is expanded in powers of the elementary symmetric functions 
E 2 , E 3 , E 4 , E 5 of X, Y, Z, P, P (since E x = ). 

ALGORITHM FOR Rj . Let x, y, z have nonnegative real part and at most one 
of them be 0, while Rep > 0. Alternatively, if p ^ and |php| < n, either let x,y,z 
be real and nonnegative and at most one of them be 0, or else let two of the variables 
x, y, z be nonzero and conjugate complex with phase less in magnitude than n and the 
third variable be real and nonnegative. The function Rj(x,y, z,p) defined by (j^j is to be 
computed with relative error less in magnitude than r. (We assume r < 10~ A .) 
Let (x ,y ,z ,p ) = (x,y,z,p) and 
. x + y + z + 2p 

A = »— S={p-x){p-y){p-z), 24 

5 

Q = (r/4)- 1 / 6 max{\A - x\, \A - y\, \A - z\, \A Q - p\} . (25) 
For m = 0,1,2, ... , define 

\Z*^m \/ Vrn y/%m \J \fy~m \J %>m i ^m+l T i 

•£"m _ Urn _ _ Pm ~t~ 

%m+l — ^ , Um+1 — ~ , z m+l — ~ , Pm+1 ~ ^ ; \& < ) 

dm i^y/Pm y/%m ) ( \J Pm ~\~ \J ym ) ( \J Pm ~\~ \J %m ) j &m "To ) (^®) 

where each square root has nonnegative real part. 

Compute A m for m — 0, 1, . . . , n , where 4r n Q < \A n \ . Compute also Rc(l, 1 + e m ) 
with relative error less in magnitude than r for m = 0, 1, . . . ,n — 1. Define 

X = 4^, ^ = 4^> ^ = 4^> P=(-X-Y-Z)/2, (29) 
4 n A n 4 n A n 4 n A n v " K ' 

E 2 = XY + XZ + YZ - 3P 2 , E 3 = XYZ + 2E 2 P + 4P 3 , (30) 
E 4 = (2XYZ + E 2 P + 3P 3 )P , E 5 = XYZP 2 . (31) 
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Then 

RAx, y, z,p) ss A~ n A- 3 / 2 ( 1 - — E 2 + ~E 3 + —El - —E 4 - —E 2 E 3 + —E 5 
jk ,y, ,yj n y 14 6 88 2 22 52 3 26 

n— 1 4— m 

+ 6 Yl -r- Rc{h 1 + e m ) (32) 

m=0 ^ m 

with relative error less in magnitude than r. 

When the variables are complex, the bound on relative error is not rigorous because 
of the possibility of some cancellation between terms that individually have error less 
than r. In practice, however, the error is usually much smaller than the bound. 

If x, y, z are real and nonnegative, at most one of them is 0, and the fourth variable 
of Rj is negative, the Cauchy principal value is given by Jl6], (4.6)] and (pi]) : 

(y + q)Rj(x,y,z,-q) = (p - y)Rj{x,y, z,p) - 3R F (x,y, z) 

( \ 1/2 

+ 3 R c (xz + pq,pq), q > 0, (33) 

\xz + pq J 

where p — y = (z — y)(y — x)/(y + q). If we permute the values of x,y,z so that 
(z — y)(y — x) > 0, then p > y > and all terms on the right side can be computed 
by the preceding algorithms. If x, y, z are complex, conditions of validity have not been 
established for this method of computing the Cauchy principal value. 

It would be desirable to give the algorithm less restrictive conditions on x, y, z,p that 
do not rule out so many cases in which p = z and Rj reduces to Rp, for in such cases 
the algorithm gives correct results under the much weaker conditions stated below in the 
algorithm for R D . We note that p = z implies e m = and d m = 2J~z^(z m + A m ) for all 



m, whence ( p3|) becomes 



n-l 



4~m 

Rd(x , y , Zq) = 4~ n i? D (x n , y n , z n ) + 3 ^ —j=- ( r v • (34) 

The function R D (x,y,z) is treated as a symmetric function of x } y , z } z } z. 



ALGORITHM FOR Rd . Let x,y,z be nonzero and have phase less in magnitude 
than 71, with the exception that at most one of x and y may be 0. The function Rd(x, y, z) 
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defined by (0j is to be computed with relative error less in magnitude than r. (We assume 
r < 10" 4 .; 

Let xq — x, y — y, z = z, and 

= x + y + 3z = (r/4) _i/ 6 max{ | Ao _ x \ {Ao _ | | Aq _ z{} (35) 

5 

For 771 = 0, 1, 2, ... , define 

y/^Cm ~\J Vm \J ■Em y/^m \fy~m y/^m j T j (*^) 

Xm+1 = T , 2/m+l = ^ , Z m+ i = , [6 I ) 

where each square root has nonnegative real part. 

Compute A m for m = 0, 1, . . . , n , where 4~ n Q < \A n \ . Define 

X = 4^, F=4^ ; Z=-(X + Y)/3, (38) 

E 2 = XY - 6Z 2 , E 3 = (3XY - 8Z 2 )Z , (39) 
E A = 3(XY - Z 2 )Z 2 , E 5 = XYZ 3 . (40) 

Then 

R D (x, y, z) « 4-"A- 3 / 2 ( 1 - — £ 2 + -£ 3 + — El - —E A - — E 2 E 3 + —E 5 ) 
v\ ,y, j n y 14 2 6 88 2 22 52 d 26 / 

n— 1 A_— m 

+ 3E 7^ +A ^ (41) 

m=0 V m r"> "my 

ifit/i relative error less in magnitude than r. 

The remark about relative error immediately following the algorithm for Rj applies 
again here. The function Rq can be computed from R F and Rd by using (|7D if at 
most one of x,y,z is 0. Neither ([/]) nor the next algorithm can be used to compute 
R G {0,0,z) = v/i/2. 

A faster and simpler way of computing the complete case of Rf and Rg with complex 
variables is useful because Legendre's complete integrals of the first and second kinds are 

K(k) = R F (l-k 2 , 1, 0), E[k) = 2R G (1 — k 2 , 1 , 0). (42) 
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The method of successive arithmetic and geometric means (a special case of successive 
Landen transformations) serves this purpose for complex k 2 ^ [1, +00). (Algorithms us- 
ing Landen transformations for incomplete R F and R G are given in || but only for real 
variables; conditions of validity for complex variables are not obvious.) 

ALGORITHM FOR R F (x , y, 0) AND R G {x , y, 0). Let x and y be nonzero and 
have phase less in magnitude than ix. The complete elliptic integrals R F (x , y , 0) and 
Rg(x , y , 0) are to be computed with relative error less in magnitude than r. 
Define x Q = a/x, y = ^fy , and 

•Em 1 ym , n 1 o / A n\ 

X m +1 = ~ , Vm+1 = V X mVm , 771 = 0,1,2,..., (43) 

where each square root has positive real part. Compute x m and y m for m = 0,1,. . . ,n, 
where 

\x n ~y n \< 2.7^/r\x n \ . (44) 



Then 



R F {x,y,0)n—?— (45) 

•^n ~r yn 



with relative error less in magnitude than r. Also, 

2R G (x,y,0)^ ^^±^) 2 - j22 m -\x m -y m f^R F {x,y^) (46) 

with relative error less in magnitude than r if we neglect terms of order r 2 . The summa- 
tion is empty if n = . 



This algorithm can be used to compute also 

Rd(0, y, z) = 3 _ [2RG{y, z,0)-z R F {y, z,0)), O^y^z^O. (47) 
The exceptional case with y = z ^ is 

R D (0,y,y) = 3 -^y- 3/2 . (48) 
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3 Numerical checks 



Codes based on the algorithms of Section 2 can be checked against the following as- 
sortment of numerical values for complete and incomplete integrals with real, conjugate 
complex, or nonconjugate complex variables. 



Rp 
Rp 
Rp 
Rp 
Rp 
Rp 

Rc 
Rc 
Rc 
Rc 
Rc 
Rc 

Rj 
Rj 
Rj 
Rj 
Rj 
Rj 
Rj 
Rj 



1,2,0) = 1.3110 28777 1461 

i, -i, 0) = R F (0.5, 1, 0) = 1.8540 74677 3014 

i - 1, i, 0) = 0.79612 58658 4234 - i 1.2138 56669 8365 

2,3,4) = 0.58408 28416 7715 

i,-i,2) = 1.0441 44565 4064 

i-l,i,l-i) = 0.93912 05021 8619 - i 0.53296 25201 8635 

0, 1/4) = 7T = 3.1415 92653 5898 
9/4, 2) = In 2 = 0.69314 71805 5995 
0, i) = (1 - i) 1.1107 20734 5396 

= 1.2260 84956 9072 - i 0.34471 13698 8768 



1/4,-2) 



In 2 



0.23104 90601 8665 



1 , -1) = 0.77778 59692 0447 + i 0.19832 48499 3429 

0. 1, 2, 3) = 0.77688 62377 8582 
2,3,4,5) = 0.14297 57966 7157 

2, 3, 4, -1 + i) = 0.13613 94582 7771 - i 0.38207 56162 4427 

1, -i,0,2) = 1.6490 01166 2711 

-1 + i, -1 - i, 1, 2) = 0.94148 35884 1220 

i, -i, 0, 1 - i) = 1.8260 11522 9009 + i 1.2290 66190 8643 

-1 + i, -1 - i, 1, -3 + i) = -0.61127 97081 2028 - i 1.0684 03839 0007 

-1 +i,-2- i, -i, -l + i) = 1.8249 02739 3704 - i 1.2218 47578 4827 
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The last case does not fit the assumptions, but the algorithm yields a value agreeing with 
Rd{— 2 — i, —i, below. If x, y, z are strictly positive, the Cauchy principal value 

of Rj changes sign at a negative value of the fourth variable, as illustrated by 

Rj(2, 3, 4, -0.5) = 0.24723 81970 3052 
Rj{2, 3, 4, -5) = -0.12711 23004 2964 

In this example Rj vanishes when the fourth variable is approximately —1.2552, and 
cancellation between terms on the right side of (|33"D will lead to loss of significant figures. 

R D (0, 2, 1) = 1.7972 10352 1034 
R D (2, 3, 4) = 0.16510 52729 4261 
R D (i, -i, 2) = 0.65933 85415 4220 

R D (0,i, -i) = 1.2708 19627 1910 + i 2.7811 12015 9521 
R D (0, i - 1, i) = -1.8577 23543 9239 - i 0.96193 45088 8839 
R D {-2 - i, -i, -1 +i) = 1.8249 02739 3704 - % 1.2218 47578 4827 

R G (0, 16, 16) = 2E(0) = ir = 3.1415 92653 5898 
Rq{2, 3, 4) = 1.7255 03028 0692 
R G (0,i, -i) = 0.42360 65423 9699 

Rg{% - 1, i, 0) = 0.44660 59167 7018 + i 0.70768 35235 7515 
R G {-i, i-\,i) = 0.36023 39218 4473 + i 0.40348 62340 1722 
^(0,0.0796,4) = ^(0.99) = 1.0284 75809 0288 

Consistency checks that do not use external information are furnished by the following 
equations, in which x,y,p are positive, A/i = xy, and |phA| < n : 

R F (x + A, y + A, A) + Rf(x + /x, y + /i, /x) = R F (x, y, 0) ; (49) 
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Rc(X, x + \) + R c (^x + fi)= R c {0, x) ; (50) 
Rj(x + \,y + \,\,p + X) + Rj(x + fj,,y + (j,,fj,,p + fj) = Rj(x, y, 0, p) - 3 R c (a, b) , (51) 
where 



a = p(\ + fi + x + y), b = p{p + \)ip + fi) , b — a — p(p — x)ip — y) ; (52) 

3 

R D (\, x + A, y + A) + R D (fi, x + fi, y + fi) = R D iO, x, y) = . (53) 

y + y + A + /i 

These equations are special cases of addition theorems. Another consistency check is 
furnished by 

3 

Rd (x, y, z) + R D (y, z, x) + R D (z, x, y) = , (54) 

where x, y, z lie in the cut plane and each square root has positive real part. 



4 Other integrals 

Legendre's complete elliptic integrals K and E are given by 

K(k) = R F iO,l-k 2 ,l), (55) 
E(k) = 2i?c(0, 1 — k 2 , 1) 

= ^^[R D (0,l-k 2 ,l) + R D (0,l,l-k% (56) 

K{k)-E{k) = ^Roi^l-k 2 ,!), (57) 

E(k) — (1 — k 2 )K(k) = fc2(1 ~ fc2) i? D (0,l,l-fc 2 ). (58) 

In Legendre's incomplete integrals we shall use the abbreviation c = csc 2 = l/sin 2 0: 

F{<f>,k) = ism <p) R F icos 2 <p, I -k 2 sin 2 (j),l) = R F ic- I, c-k 2 ,c), (59) 

k 2 

E{<j>,k) = R F ic-l,c-k 2 ,c)-—R D ic-l,c-k 2 ,c), (60) 



n(0,A;,n) = / (l + nsiii 2 ^)- 1 ^ - k 2 sin 2 6)- 1 ' 2 d6 
Jo 



= _Rf(c — 1, c — k 2 , c) — — Rjic — 1, c — A; 2 , c, c + n) . (61) 
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Some related integrals are 

D((f),k) = sin 2 6(l-k 2 sin 2 6)- 1/2 d6 = - R D (c - 1, c - k 2 , c) , (62) 
Jo 3 

k 2 i 

K(k)Z(P,k) = — sm(3cos(3^Jl-k 2 sin 2 f3Rj{0,l-k 2 ,l,l-k 2 sm 2 f3), (63) 
3 



A ((3,k) 



2 (1 - k 2 ) sin $ cos /3 



7T 



A 



k 2 ( k 2 \ 
R F (0, 1 - k 2 , 1) + — Rj f 0, 1 - /c 2 , 1, 1 - —J 



(64) 



where A = y 1 — (1 — k 2 ) sin 2 /3. The function Z((3, k) is Jacobi's zeta function p, 140.03], 
and A (/3,k) is Heuman's lambda function p, 150.01]. 
Bulirsch's elliptic integrals |ij are 



el l(x, k c ) 
el 2(x, k c , a, b) 



= x R^l.l + k 2 c x 2 ,1 + x 2 ) , 
= ax R F {l,l + k 2 c x 2 ,1 + x 2 ) 
+ hb-a)x 3 R D (l,l + k 2 c x 2 ,l + x 2 ) 



el3(x,k C} p) = xRp(l, 1 + k c x , 1 + x 



cel(k c ,p, a, b) 



+ -(1 - p)x 3 Rj(l, 1 + k 2 c x 2 , 1 + x 2 , 1 + px 2 ) 
= a R F (0, k 2 c , 1) + -(6 - pa)Rj(0, k 2 , l,p) . 



(65) 

(66) 

(67) 
(68) 



In the real domain, inverse circular and inverse hyperbolic functions are expressed in 
terms of Rc by 

' x\ 



m'~J = (a; - y)Rc ((~y^) 



arctan(a;/|/) = xRciy 2 ,y 2 + x 2 ), 

arctanh(x/y) = xRc(y 2 ,y 2 — x 2 ), 

&rcsm(x/y) = xRciy 2 — x 2 ,y 2 ), 

arcsinh(x/?/) = xRc(y 2 + x 2 ,y 2 ), 

arccos(x/y) = (y 2 - x 2 )^R c (x 2 , y 2 ). 

arccosh(x/?/) = (x 2 — y 2 )^R c (x 2 ,y 2 ), 



x > 0, 

— oo < x < oo, 
-y < x < y, 
-y <x <y, 

— OO < X < oo, 

< x < y, 
x>y, 



(69) 
(70) 
(71) 
(72) 
(73) 
(74) 
(75) 



where y > in each case. These equations remain valid in the complex domain provided 
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that the variables of Rc satisfy the conditions accompanying (|3]). If y = 1 the function 
multiplying Rq shows in each case the asymptotic behavior as the left side tends to 0. 
Many elliptic integrals of the form 

rf[{ai + bit)*/ 2 dt, (76) 
J y i=i 



where pi,...,p n are integers and the integrand is real, are reduced in |]7|-pT| to the 
integrals in Section 2. Use of the algorithm for Rp is illustrated by numerical computation 
of the integral 

i=r , dt (77) 

J « \]Ui + 2ji« + M 2 )(/ 2 + 2g 2 t + h 2 e) 
where all quantities are real, x > y, the two quadratic (or, if hi = 0, linear) polynomials 
are positive on the open interval of integration, and their product has at most simple 
zeros on the closed interval. Let 

qi {t) = f i + 2g i t + h i t 2 , i = l,2, (78) 



(x-y)U = yjqi(x)<b(y) + ^q 1 (y)q 2 (x) , (79) 
T = 2g 1 g 2 -hh 2 -foh 1 , (80) 



V = 2J(g!-f 1 h 1 )(g 2 2 -f 2 h 2 ). (81) 



Then 



1 = 2 R F (U 2 + T + V, U 2 + T -V, U 2 ) . (82) 
Except for notation this is the same as @, (34)]. If the interval of integration is infinite, U 



is obtained by taking a limit; for example, if x = +oo , then U = \Jhiq 2 {y) + \J qi(y)h 2 ■ If 
exactly one of q\ and q 2 has conjugate complex zeros, then V is pure imaginary; otherwise 
the variables of Rp are real. The algorithm for Rp in this paper can be used in both 
cases. However, if both polynomials have conjugate complex zeros, the quadrilateral 
with the zeros as vertices has diagonals that must not intersect at an interior point of 
the interval of integration. If they do, the integral must be split into two parts at the 
point of intersection. This restriction is discussed in M, §4] and removed by a Landen 



transformation in 11 . 
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